In RodolfoPelinson/pelinson.et.al.2020: Data Analysis of Pelinson et al 2020

knitr::opts_chunk$set(echo = TRUE)

The goal of pelinson.et.al.2020 is to walk the user through the statistical analysis presented in:
"Pelinson et al 2020. Top predator introduction changes the effects of spatial isolation on freshwater community structure"
DOI: https://doi.org/10.1101/857318

You can install the last version of PredatorIsolationComm package from my GitHub with:

``` {r installing package not eval, eval = FALSE, echo = T} install.packages("devtools") devtools::install_github("RodolfoPelinson/PredatorIsolationComm") library("PredatorIsolationComm")

This will give you access to all the data and functions used to produce the results shown in "Pelinson et al 2020. Top predator introduction changes the effects of spatial isolation on freshwater community structure".


``` {r installing package, eval = T, echo = FALSE}
library("PredatorIsolationComm")

Other packages used here are:
lme4 version 1.1-23
emmeans version 1.4.8

``` {r installing packages2, eval = T, results = "hide", warning = F, message = F} library(lme4) library(emmeans)

## Testing for differences in abundance

First, lets load the necessary data.
```r
data(abundance_predators)
data(abundance_consumers)
data(survey)
data(fish)
data(isolation)
data(ID)

This is for testing for differences in the abundance of predatory insects across all treatments and sampling surveys. We used generalized linear mixed models with a negative binomial distribution (glmer.nb function from package lme4) to fit the models.

First we looked for the best probability distribution to model our abundance data. We considered Gaussian, Poisson and Negative Binomial distributions.

NB_model <- glmer.nb(abundance_predators~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))

P_model <- glmer(abundance_predators~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"), family = "poisson")

G_model <- lmer(abundance_predators~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = lmerControl(optimizer = "bobyqa"), REML = F)

plot(G_model, main = "Gaussian")
plot(P_model, main = "Poisson")
plot(NB_model, main = "Negative Binomial")

AIC(G_model)
AIC(P_model)
AIC(NB_model)

We chose to use the Negative Binomial distribution to model data after inspecting AIC values (Negative Binomial had the lowest value) for each model and the spread of Pearson residuals against fitted values.

Then we used anova from package lme4 to compute likelihood ratio tests.

`pred_no_effect` <- glmer.nb(abundance_predators~ 1 + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_survey` <- glmer.nb(abundance_predators~(survey) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_fish` <- glmer.nb(abundance_predators~(survey+fish) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_isolation` <- glmer.nb(abundance_predators~(survey+fish+isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_survey:fish` <- glmer.nb(abundance_predators~(survey + fish + isolation + survey:fish) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_survey:isolation` <- glmer.nb(abundance_predators~(survey + fish + isolation + survey:fish + survey:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_fish:isolation` <- glmer.nb(abundance_predators~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`pred_survey:fish:isolation` <- glmer.nb(abundance_predators~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))

Anova_predators <- anova(`pred_no_effect`, 
      `pred_survey`,
      `pred_fish`,
      `pred_isolation`,
      `pred_survey:fish`,
      `pred_survey:isolation`,
      `pred_fish:isolation`,
      `pred_survey:fish:isolation`)

Anova_predators

Now some post-hoc tests using function emmeans from package emmeans to identify pairwise differences for the effect of sampling survey and isolation, separately. Post-hoc tests were always applied to the most complex model to account for the effect of all possible interactions.

emmeans(`pred_survey:fish:isolation`, list(pairwise ~ survey), adjust = "sidak") 
emmeans(`pred_survey:fish:isolation`, list(pairwise ~ isolation), adjust = "sidak")

Same thing now for herbivores and detritivores.
Again, we first looked for the best probability distribution to model our abundance data. We considered Gaussian, Poisson and Negative Binomial distributions.

NB_model <- glmer.nb(abundance_consumers~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))

P_model <- glmer(abundance_consumers~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"), family = "poisson")

G_model <- lmer(abundance_consumers~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = lmerControl(optimizer = "bobyqa"), REML = F)

plot(G_model, main = "Gaussian")
plot(P_model, main = "Poisson")
plot(NB_model, main = "Negative Binomial")

AIC(G_model)
AIC(P_model)
AIC(NB_model)

We, again, chose to use the Negative Binomial distribution to model data after inspecting AIC values (Negative Binomial had the lowest value) for each model and the spread of Pearson residuals against fitted values.

`cons_no_effect` <- glmer.nb(abundance_consumers~ 1 + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_survey` <- glmer.nb(abundance_consumers~(survey) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_fish` <- glmer.nb(abundance_consumers~(survey+fish) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_isolation` <- glmer.nb(abundance_consumers~(survey+fish+isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_survey:fish` <- glmer.nb(abundance_consumers~(survey + fish + isolation + survey:fish) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_survey:isolation` <- glmer.nb(abundance_consumers~(survey + fish + isolation + survey:fish + survey:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_fish:isolation` <- glmer.nb(abundance_consumers~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))
`cons_survey:fish:isolation` <- glmer.nb(abundance_consumers~(survey + fish + isolation + survey:fish + survey:isolation + fish:isolation + survey:fish:isolation) + (1|ID), data = treatments, control = glmerControl(optimizer = "bobyqa"))

Anova_consumers <- anova(`cons_no_effect`, 
      `cons_survey`,
      `cons_fish`,
      `cons_isolation`,
      `cons_survey:fish`,
      `cons_survey:isolation`,
      `cons_fish:isolation`,
      `cons_survey:fish:isolation`)
Anova_consumers

Now post-hoc tests to identify pairwise differences for the effect of sampling survey, and interaction between isolation and survey, and between isolation and presence of fish:

emmeans(`cons_survey:fish:isolation`, list(pairwise ~ survey), adjust = "sidak") 
emmeans(`cons_survey:fish:isolation`, list(pairwise ~ isolation|fish), adjust = "sidak") 
emmeans(`cons_survey:fish:isolation`, list(pairwise ~ isolation|survey), adjust = "sidak")

Finaly, we can plot the abundances of predators and herbivores/detritivores for each survey.

abundance_SS1_predators <- abundance_predators[which(survey == "1")]
abundance_SS2_predators <- abundance_predators[which(survey == "2")]
abundance_SS3_predators <- abundance_predators[which(survey == "3")]

abundance_SS1_consumers <- abundance_consumers[which(survey == "1")]
abundance_SS2_consumers <- abundance_consumers[which(survey == "2")]
abundance_SS3_consumers <- abundance_consumers[which(survey == "3")]

par(mfrow = c(1,2))

boxplot(abundance_SS1_predators~fish_isolation_SS1, outline = F, ylab = "Abundance", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, ylim = c(0,240), col = "transparent", main = "First Survey - Predators", xaxt="n")
mylevels <- levels(fish_isolation_SS1)
levelProportions <- summary(fish_isolation_SS1)/length(fish_isolation_SS1)
col <- c(c("sienna1","sienna3","sienna4"), c("steelblue1","steelblue3","steelblue4"))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){

  x<- c(1,2,3,5,6,7)[i]
  thislevel <- mylevels[i]
  thisvalues <- abundance_SS1_predators[fish_isolation_SS1==thislevel]

  # take the x-axis indices and add a jitter, proportional to the N in each level
  myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.8)
  points(myjitter, thisvalues, pch=pch[i], col=col[i] , cex = 1.5, lwd = 3) 

}
boxplot(abundance_SS1_predators~fish_isolation_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n")
axis(1,labels = c("30m", "120m", "480m","30m", "120m", "480m"), cex.axis = 0.8, at =c(1,2,3,5,6,7))
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
box(lwd = 2.5)


boxplot(abundance_SS1_consumers~fish_isolation_SS1, outline = F, ylab = "Abundance", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, ylim = c(0,1250), col = "transparent", main = "First Survey - Non Predators", xaxt="n")
mylevels <- levels(fish_isolation_SS1)
levelProportions <- summary(fish_isolation_SS1)/length(fish_isolation_SS1)
col <- c(c("sienna1","sienna3","sienna4"), c("steelblue1","steelblue3","steelblue4"))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){

  x<- c(1,2,3,5,6,7)[i]
  thislevel <- mylevels[i]
  thisvalues <- abundance_SS1_consumers[fish_isolation_SS1==thislevel]

  # take the x-axis indices and add a jitter, proportional to the N in each level
  myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.8)
  points(myjitter, thisvalues, pch=pch[i], col=col[i] , cex = 1.5, lwd = 3) 

}
boxplot(abundance_SS1_consumers~fish_isolation_SS1, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n")
axis(1,labels = c("30m", "120m", "480m","30m", "120m", "480m"), cex.axis = 0.8, at =c(1,2,3,5,6,7))
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
box(lwd = 2.5)


boxplot(abundance_SS2_predators~fish_isolation_SS2, outline = F, ylab = "Abundance", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, ylim = c(0,240), col = "transparent", main = "Second Survey - Predators", xaxt="n")
mylevels <- levels(fish_isolation_SS2)
levelProportions <- summary(fish_isolation_SS2)/length(fish_isolation_SS2)
col <- c(c("sienna1","sienna3","sienna4"), c("steelblue1","steelblue3","steelblue4"))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){

  x<- c(1,2,3,5,6,7)[i]
  thislevel <- mylevels[i]
  thisvalues <- abundance_SS2_predators[fish_isolation_SS2==thislevel]

  # take the x-axis indices and add a jitter, proportional to the N in each level
  myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.8)
  points(myjitter, thisvalues, pch=pch[i], col=col[i] , cex = 1.5, lwd = 3) 

}
boxplot(abundance_SS2_predators~fish_isolation_SS2, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n")
axis(1,labels = c("30m", "120m", "480m","30m", "120m", "480m"), cex.axis = 0.8, at =c(1,2,3,5,6,7))
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
box(lwd = 2.5)



boxplot(abundance_SS2_consumers~fish_isolation_SS2, outline = F, ylab = "Abundance", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, ylim = c(0,1250), col = "transparent", main = "Second Survey - Non Predators", xaxt="n")
mylevels <- levels(fish_isolation_SS2)
levelProportions <- summary(fish_isolation_SS2)/length(fish_isolation_SS2)
col <- c(c("sienna1","sienna3","sienna4"), c("steelblue1","steelblue3","steelblue4"))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){

  x<- c(1,2,3,5,6,7)[i]
  thislevel <- mylevels[i]
  thisvalues <- abundance_SS2_consumers[fish_isolation_SS2==thislevel]

  # take the x-axis indices and add a jitter, proportional to the N in each level
  myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.8)
  points(myjitter, thisvalues, pch=pch[i], col=col[i] , cex = 1.5, lwd = 3) 

}
boxplot(abundance_SS2_consumers~fish_isolation_SS2, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n")
axis(1,labels = c("30m", "120m", "480m","30m", "120m", "480m"), cex.axis = 0.8, at =c(1,2,3,5,6,7))
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
box(lwd = 2.5)



boxplot(abundance_SS3_predators~fish_isolation_SS3, outline = F, ylab = "Abundance", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, ylim = c(0,240), col = "transparent", main = "Third Survey - Predators", xaxt="n")
mylevels <- levels(fish_isolation_SS3)
levelProportions <- summary(fish_isolation_SS3)/length(fish_isolation_SS3)
col <- c(c("sienna1","sienna3","sienna4"), c("steelblue1","steelblue3","steelblue4"))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){

  x<- c(1,2,3,5,6,7)[i]
  thislevel <- mylevels[i]
  thisvalues <- abundance_SS3_predators[fish_isolation_SS3==thislevel]

  # take the x-axis indices and add a jitter, proportional to the N in each level
  myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.8)
  points(myjitter, thisvalues, pch=pch[i], col=col[i] , cex = 1.5, lwd = 3) 

}
boxplot(abundance_SS3_predators~fish_isolation_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n")
axis(1,labels = c("30m", "120m", "480m","30m", "120m", "480m"), cex.axis = 0.8, at =c(1,2,3,5,6,7))
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
box(lwd = 2.5)



boxplot(abundance_SS3_consumers~fish_isolation_SS3, outline = F, ylab = "Abundance", xlab = "", at = c(1,2,3,5,6,7), lwd = 1.5, ylim = c(0,1250), col = "transparent", main = "Third Survey - Non Predators", xaxt="n")
mylevels <- levels(fish_isolation_SS3)
levelProportions <- summary(fish_isolation_SS3)/length(fish_isolation_SS3)
col <- c(c("sienna1","sienna3","sienna4"), c("steelblue1","steelblue3","steelblue4"))
#bg <- c(rep("sienna3",3), rep("dodgerblue3",3),rep("sienna3",3), rep("dodgerblue3",3))
pch <- c(15,16,17,15,16,17)
for(i in 1:length(mylevels)){

  x<- c(1,2,3,5,6,7)[i]
  thislevel <- mylevels[i]
  thisvalues <- abundance_SS3_consumers[fish_isolation_SS3==thislevel]

  # take the x-axis indices and add a jitter, proportional to the N in each level
  myjitter <- jitter(rep(x, length(thisvalues)), amount=levelProportions[i]/0.8)
  points(myjitter, thisvalues, pch=pch[i], col=col[i] , cex = 1.5, lwd = 3) 

}
boxplot(abundance_SS3_consumers~fish_isolation_SS3, add = T, col = "transparent", outline = F,at = c(1,2,3,5,6,7), lwd = 1.5, xaxt="n")
axis(1,labels = c("30m", "120m", "480m","30m", "120m", "480m"), cex.axis = 0.8, at =c(1,2,3,5,6,7))
axis(1,labels = c("Fishless","Fish"), cex.axis = 1, at =c(2,6), line = 1.5, tick = F )
box(lwd = 2.5)